Quantum quench influenced by an excited-state phase transition 
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We analyze excited-state quantum phase transitions (ESQPTs) in three schematic (integrable and 
nonintegrable) models describing a single-mode bosonic field coupled to a collection of atoms. It is 
shown that the presence of the ESQPT in these models affects the quantum relaxation processes 
following an abrupt quench in the control parameter. Clear cut evidence of the ESQPT effects is 
presented in integrable models, while in the nonintegrable model the evidence is blurred due to 
chaotic behavior of the system in the region around the critical energy. 
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I. INTRODUCTION 

Diverse quantum effects in systems depending on ex- 
ternal control parameters represent an interesting field of 
theoretical and experimental investigation. A lot of re- 
cent attention in this field has been focused on two differ- 
ent types of dynamical phenomena, namely on quantum 
phase transitions and so-called quantum quenches. These 
phenomena and their mutual relation are addressed in 
the present work. 

A quantum phase transition (QPT) is a sudden change 
of the ground-state structure at a certain critical value 
of the control parameter A. ft can be observed as a non- 
analytic evolution of the system's energy and wave func- 
tion induced by an adiabatic variation of the control pa- 
rameter A across the quantum critical point A c o at zero 
temperature. First discussed in the 1970s the QPT 

phenomena become very important in the context of solid 
state physics as well as in nuclear and many-body 
physics — see, e.g., recent reviews @, 

A quantum quench (QQ) represents an abrupt, dia- 
batic change Ai — > A2 of the control parameter followed 
by a system-specific quantum relaxation process. Pio- 
neering theoretical works in this field appeared already 
in the late 1960s @, but a really rapid growth of interest 
was triggered by experimental studies at the beginning of 
this millennium [lOl. For an extensive list of QQ-related 
references see Ref. [111 ]. 

Not surprisingly, the QPT and QQ effects can be mu- 
tually related. If the initial state before a quench coin- 
cides with the ground state near a QPT, the dynamics 
after the quench depends substantially on whether the 
parameter change does or does not bring the system to 
the other quantum phase, or eventually to the narrow 
quantum critical region between the phases [Ill - ll4j ]. 

In this paper, we discuss the quench-induced dynamics 
in connection with a novel concept related to quantum 



criticality — a so called excited-state quantum phase tran- 
sition (ESQPT) [l5l - [l7j |. This phenomenon represents a 
nonanalytic evolution of individual excited states in the 
system with a variable control parameter. So far, such 
effects have been studied mostly in integrable systems 
with one effective degree of freedom (one dimensional 
configuration spaces), showin g, a sing ularity of classical 
dynamics at a certain energy [18l425l |. but they seem to 
exist in a much richer variety of incarnations. 

The ESQPTs can be viewed as a reinterpretation of 
thermal phase transitions in the microcanonical lan- 
guage. The key step is a scaling of energy and other 
observables by a suitably defined size parameter H. Let 
(•)t stands for a thermal average of the quantity in 
brackets at temperature T. If the thermal fluctuation 
(AS 2 ) T = (£ 2 )t - (£)t of the scaled energy £ = E/X 
vanishes in the thermodynamic limit H — > 00 (often syn- 
onymous with the classical limit h — > 0) , a thermal phase 
transition at temperature T = T c becomes localized at a 
sharp value £ c = (£)t c of the scaled energy. This shows 
up as an anomalous (nonanalytic) "flow" of energy levels 
£i(X) through the boundary £ C (X) defining the locus of 
ESQPT points [H-Ql H . 

Besides the shape of individual £i(X) curves, the 
excited-state transition affects also the dependence of the 
level density p(£) on the scaled energy. The function p(£) 
at £ — £ c shows a singularity whose type depends on the 
underlying thermal phase transition and/or on the corre- 
spondin g an omaly in the associated classical phase space 
[3 fl8L I19I] . Since the singularity influences quantum 
relaxation processes in the critical region, we can antici- 
pate a major impact of the ESQPT on the QQ-induccd 
dynamics in those cases in which the energy distribution 
after the quench is centered near £ c . 

To study the above-formulated conjecture, we use three 
simple models. They describe a single-mode bosonic field 
interacting with an algebraic subsystem, which is based 
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either on the SU(1,1) or on the SU(2) dynamical algebra. 
The SU(1,1) model may serve as a toy for the description 
of formation and dissociation of diatomic molecules and 
bosonic atoms [26J. The SU(2) Hamiltonian represents 
either the well-known Dicke [27| or Jaynes-Cummings 
[iH (Tavis-Cummings (29J) models of quantum optics, 
or may alternatively describe an interacting mixture of 
diatomic molecules and fermionic atoms [261 ] . In this 
work we will adopt the former interpretation. While the 
SU(1,1) model and the SU(2)-based Jaynes-Cummings 
model are integrable, the SU(2)-based Dicke model is not. 

All the three models show a rather similar phase struc- 
ture. First, considering the ground-state properties, 
it turns out that an increasing strength of interaction 
between the bosonic field and the algebraic subsystem 
drives the entire system to a quantum critical point where 
the ground state abruptly changes its form. We then 
show that this QPT is in all three models followed by 
a chain of ESQPTs and demonstrate that these have a 
strong impact on the character of relaxation dynamics 
after some fine-tuned quantum quenches — namely those 
leading the system to a narrow region around the critical 
excitation energy. The type of the ESQPT and its QQ 
signatures depends on the dimensionality of the model: 
they are of the strongest type (see below for a precise 
definition of this concept) for the SU(1,1) and the SU(2) 
Jaynes-Cummings Hamiltonians, and of a softer type for 
the SU(2) Dicke Hamiltonian. We will also see that 
breaking of integrability in the latter model blurs the ef- 
fects of criticality in the quench dynamics. We anticipate 
the same general trend also in more complex situations. 

The plan of the paper is the following: In Sec. [XT] we 
describe the models, and analyze in Sec. Mil their classi- 
cal and phase-transitional properties, particularly those 
related to excited states. Sec. IIVI collects the results on 
quantum quenches. We introduce a general concept of a 
critical quench, driving the system into the ESQPT re- 
gion, and continue to more specific numerical results for 
the models employed. Sec.fVlbrings a brief summary and 
outlook. 



II. MODELS 

A. Algebraic structure 

Below we will investigate quantum quenches in two 
simple models, both describing a system composed of 
two interacting parts: (i) a single bosonic mode given by 
creation and annihilation operators and b, therefore 
described by the Heisenberg-Weyl algebra HW(1), and 
(ii) a system represented by pseudospin operators J± — 
J x ±iJ y and Jo = J z satisfying commutation relations of 
the SU(2) algebra, 

[J 0) J±]=±J±, [J +) J_] = 2J 0) (1) 




FIG. 1: A schematic representation of the models used. The 
SU(1,1) model (the upper panel) may describe the coexistence 
of two-atom molecules (lower level) with dissociated atoms 
(upper level). The SU(2) model (the lower panel) describes 
the interaction of a single-mode radiation field with an array 
of two-level atoms. 

or by analogous operators K± = K x ± iK y and Kq = K z 
satisfying commutation relations of the SU(1,1) algebra 

[K , K±] = ±K±, [K+, AL] = -2K . (2) 

The SU(2) or SU(1,1) algebras will be realized more 
specifically in terms of fermionic or bosonic operators. 
The complete dynamical algebra is HW(1) <g> SU(»), 
where the bullet stands for a specification of the respec- 
tive special unitary algebra, but in the following we will 
use just abbreviated names SU(2) and SU(1,1) for the 
two models. A schematic representation of both models 
is given in Fig. [TJ 

The Hilbert space of the coupled system is identified 
with the tensor product H = HW®H^, where is the 
space of HW(1) spanned by the set of basis vectors \Nb), 
with Nf, — 0, 1, . . . denoting the number of 6-bosons, and 
i'"' coincides with the space associated with one of the 
irreducible representations (irreps) of the groups SU(2) 
or SU(1,1). The irreps are classified by the eigenvalues 

(2) 

c su(») °f ^ ne res P ec tive second-order Casimir invariant, 

r (2) - T 2 4 T 2 4 T 2 

SU(2) — x i :J y i °z i \°) 

4u(i4) = K 2 +K 2 ~Kl (4) 

(2) 

which are parametrized as c gl j( 2 ) = j(j + 1) (with j in- 

(2) 

teger or half-integer) and Cg-yn ^ = — k(k — 1) (with 
k > known as the Bergmann index). The irreps are 
finite-dimensional in the SU(2) case (compact group) 
and infinite-dimensional in the SU(1,1) case (noncom- 
pact group). The respective basis states \j, to) (with 
m = — j, + being the eigenvalue of Jo) and 

\k,n) (with n = 0,1,... enumerating the eigenvalues 
k + n of Kq) are generated from the lowest state \j, —j) 
and |fe,0) by consecutive actions of the rising operators 
J + and K + , respectively. 
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B. SU(1,1) model 

Since the group generated by the SU(1,1) algebra is 
noncompact, its irreps are infinite dimensional, the gener- 
ators being expressible through creation and annihilation 
operators a\a of another type of bosons. For instance, 
a single boson pair realization reads as 



K+ = 5(a f ) 2 , K. 



±a 2 , K a = \ (a t a+ ±) . (5) 



Alternatively, one can use some other boson pair real- 
izations (e.g., with two kinds of bosons), which together 
with Eq. ([5]) constitute various forms of the Schwinger 
representation of the SU(1,1) algebra. Note that in 
Sec. IIII Al we will also introduce the Holstein-Primakoff 
bosonic representation. 

To construct the Hamiltonian, we assume the simplest 
realization ([5]). In this case, there are just two irreps, 
one with k = \ and the other with k = 4. Their respec- 
tive Hilbert spaces are spanned by vectors |iV a ) contain- 
ing even and odd numbers of a-bosons. The interaction 
between a- and 6-bosons is considered such that the cre- 
ation of one o-boson leads to the destruction of a pair of 
a-bosons and vice versa. Such a model can schematically 
describe, for example, the formation and dissociation of 
two-atom molecules 26] . The total Hamiltonian reads as 



H 



(i) - 



A 



bK+ + tfK. 



(6) 



where A/VmW > is a scaled coupling parameter (the 
meaning of MW will be explained below) and uj, uiq stand 
for single-particle energies (we set H = 1). 

For each of the SU(1,1) irreps (classified by the quan- 
tum number fc), there are two commuting operators 
(quantum degrees of freedom) which determine the basis 
in the whole Hilbert space of physical states: one is asso- 
ciated with the number of 6-bosons, N b = b^b, the other 
with Kq, or equivalently with the number of a-bosons 
N a = a^a. At the same time, there exist two different 
integrals of motions: one is the energy H and the other 
one can be written in the form 



M« = 2N b + N a 



4fc-l 



2{N b + K -k). (7) 



The value of M (1) > is always even, MW/2 counting 
the number of 6-bosons plus the number of a-boson pairs. 
The conservation of M^ 1 ' implies that the Hamiltonian 
(j6]) represents an integrable system, which for each fixed 
value of AfW can be associated with an effective one- 
dimensional configuration space (one quantum degree of 
freedom) . 

In the following we will assume luq > lu for the SU(1,1) 
model. This means that the A = ground state can be 
identified with a molecular condensate with no pair of 
atoms. It has the form \N b = A/W /2) ® \k, 0), where the 
first term represents a state with a maximal number of b 
bosons and the second one stands for the lowest weight 
SU(1,1) state with the minimal value of N a — (N a = 1) 



for k = \ (k = |). However, for sufficiently large values 
of the coupling parameter A the interaction between the 
molecules and atomic pairs supports a more balanced 
distribution of the expectation values (N a ) and (2N b ). 
With an increasing size of the system the crossover be- 
tween the two types of the ground-state structure is get- 
ting sharper and in the infinite-size limit, M^- 1 ' — » oo, it 
becomes a phase transition. The calculation of the criti- 
cal value of the interaction strength A^ will be presented 
in Sec. IhTbI 



C. SU(2) model 

The SU(2) algebra yields a compact group with finite- 
dimensional irreps. Its generators can therefore be con- 
structed from fermionic operators. For instance, they 
can be associated with an array of spin-i particles (or 
two-level atoms) located on 2j sites: 

2j 2j 

i=l i=l 
2j 



J o - \ (°t 

i=l 



(8) 



Here, aj^ or a^i and or create or annihilate spin- 
up and spin-down states of the fermion on site i and 
the ladder operators J± describe spin flips along the ar- 
ray. Alternatively, one can use fermion pair realizations 
of the SU(2) algebra (with J± creating and annihilat- 
ing a pair of fermions), or the Schwinger or Holstein- 
Primakoff types of bosonic realizations with truncated 
Hilbert spaces (the latter bosonic realization will be dis- 
cussed in Sec. IIII Ap . Depending on the specific real- 
ization, the model can receive different physical inter- 
pretations. Below we will implicitly consider the realiza- 
tion ([5]), which may schematically describe interactions of 
single-frequency photons with two-level atoms in maser- 
like systems. 

The Hamiltonian is taken in either of the following 
forms, 



H (2) = uj J +ujtfb + 
#( 3 ) = u Q Jo + utfb + 



A 



y/Wn L 

A 



bJ+ + & f J_ 



(6 + &t)(J_ + J+) 



(9) 



,(10) 



where A/VM® or A/vM^ is a properly scaled coupling 
parameter (A > 0) and ui,uio two single-particle energies. 
The Hamiltonian is known as the Jaynes-Cummings 
[28j or Tavis-Cummings model [29j], while the Hamilto- 
nian H ( 3 ) is referred to as the Dicke model [27] . Note that 
a so-called rotating-wave approximation of H ^ leads to 
the simpler Hamiltonian H^K 

The Jaynes-Cummings Hamiltonian (j9|) is very similar 
to that of Eq. ©. It conserves the quantity 



M< 2 ) =2(N b + J +j), 



(11) 
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analogous to Eq. (JT)), and therefore corresponds to 
an integrable system described effectively by a one- 
dimensional configuration space. (The full model has 
again two degrees of freedom, associated with commut- 
ing operators Nb and Jo.) The Dicke model violates the 
conservation of M^ 2 \ but it still conserves the parity 
IT = (— 1) M< ] I 2 labeling individual eigenstates. In this 
case, the size parameter is taken as = 4j, which is 

the total number of fermionic states (twice the number 
of sites). 

As in the SU(1,1) case, the ground states of both the 
SU(2) Hamiltonians change their nature suddenly as the 
coupling strength A increases above a certain value, the 
transition having a critical character in the infinite-size 
limit, Af( 2 ),M( 3 ) — > oo. For the Hamiltonian from 
Eq. (jnj), we will assume uj > ujq, identifying the A = 
ground state with a photon vacuum, Nb — 0, combined 
with a maximally excited state of the atom array: Jq — 
iM< 2 ' - j (below we set = 4j so that J = +j at 

A = 0). At the critical coupling strength, this structure 
eventually changes into a state with (Nb) > 0, in which 
a part of energy is transferred from atoms to the photon 
field. It should be stressed that here we are running the 
model in a nonstandard regime, taking into account only 
a finite set of states with a single fixed value of M^ 2 \ 
This is in contrast to the rotating-wave approximation 
of the Dicke model, for which one usually considers the 
infinite spectrum with all values of M^ 2 \ 

For the Dicke Hamiltonian in Eq. (fl~0|) . we set 

ujq = uj, which corresponds to the resonance absorption 
and emission of photons by the atoms. The A = ground 
state has the form \Nb = 0) ® |j, —j), describing a pho- 
ton vacuum and an unexcited array of atoms (recall that 
for the Dicke model Nb + Jo is not conserved). For a 
sufficiently strong interaction between matter and light, 
the ground state flips to a form with (Jo) > —j and 
(Nb) > 0, showing a macroscopic excitation of both sub- 
systems [3(3, HH ■ This may be considered as a toy exam- 
ple of the maser phase transition. 



D. Numerical solution 

The model Hamiltonians from Eqs. ©, © and (flU)) 
can be diagonalized numerically in an appropriate basis. 
Due to the tensor product structure of the Hilbert space, 
the basis is naturally chosen in the form \Nb) (g> \k, n) for 
the SU(1,1) model, and \N b )<2)\j, m) for both SU(2)-based 
models, where \Nb) £ stands for a state with a given 
number of b bosons while \k,n), \j,m) € B^ 11 ) are basis 
vectors of the respective SU(1,1) or SU(2) irreps. Recall 
that in the SU(1,1) case, the link of \k, n) with the states 
\N a ), counting the number of a bosons, is achieved via 
setting fc = i or | for the even- or odd-iV a irreps, re- 
spectively, and n = ^(N a — 4fc ~ 1 ). Matrix elements of 





SU(1,1) 


SU(2) 




n = 1 


n = 2 


n = 3 




N a even N a odd 


Jayn.-Cumm. 


Dicke 


UJ,UJ() 


LU0 — UJ = 1 = UJ 


u) — luq = 1 = i^o 


U) = 1 = UJo 


M (n) 


2N b +N a 2N b +N a -l 


4j 


4j 




1 3 


l 

2 


l 

2 


\(n) 
A c() 


1 

V2 


1 

%/2 


1 
V2 


£i n) 


1 
2 


1 

4 


1 

4 



TABLE I: Summary of parameter setting for the three mod- 
els employed: frequencies o;,cjo, the size parameter M (n) , a 
parameter i?*-™' in the potential energy, the QPT critical pa- 



rameter A^q\ and the ESQPT critical scaled energy £, 



(n) 



. 2V a . 2 

individual Hamiltonian terms in these bases can be easily 
calculated from the known action of the and b opera- 
tors on the vectors \Nb) and the action of {K + , Kq} 



or {J+, J_, Jo} operators on vectors \j,m) or \k,n). 

For both SU(1,1) and SU(2) integrable models, the ba- 
sis includes a finite set of vectors. These are determined 
by the chosen values of the size parameters or M^ 2 \ 
which permit only a finite number of combinations of 
(Nb,n) or (Nb 7 m) satisfying Eqs. ((TJ and (fTTj) . respec- 
tively. Thus the corresponding Hamiltonian matrices are 
finite and the diagonalization is just a routine problem. 

The nonintcgrablc Dicke model, on the other hand, has 
no conservation-dictated constraint on the allowed com- 
binations of basis vectors. Its basis is therefore infinite 
and must be numerically truncated for Nb > iVtrunc, mak- 
ing the convergence tests of the diagonalization outputs 
an important issue. Speaking qualitatively, the trunca- 
tion with a fixed iVtrunc can be safely done in the low- 
energy part of the spectrum only for a sufficiently small 
interaction strength A between the b bosons (photons) 
and atoms. Indeed, we know that for A = 0, the Dicke 
Hamiltonian is diagonal in the \Nb) ® basis, the 

states with increasing Nb being associated with increas- 
ingly high excitations. Therefore, for moderate values of 
A the states with high photon numbers are only weakly 
admixed to the low-lying energy eigenstates of the sys- 
tem. As A increases, the cutoff parameter iVtrunc must 
increase accordingly for a given low-energy portion of the 
spectrum to be well reproduced. We stress that the re- 
sults presented below were tested for stability against the 
change of iV trunc (see Sec. IIVD]) . 

The setting of the model parameters, as used in the 
calculations below, and some important model-specific 
values are summarized in Tab. U (some of the symbols 
will be explained later). 



III. PHASE TRANSITIONS 
A. Classical limit 

The above Hamiltonians are specimens of a rather 
large general class of systems — namely those described 
by finite algebraic models [32j]. For these models, the 
relevant observables are constructed in terms of a fi- 
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nite set of generators d closing a dynamical algebra 
[Gi,Gj] = J2k c ijkGk with structure constants Cy/.. The 
corresponding systems have a finite number of degrees of 
freedom and their thermodynamic (infinite size) limit co- 
incides with the classical limit h — > 0. To see this, recall 
that the thermodynamic limit is generally achieved for 
asymptotic values of a properly defined size parameter H 
such that thermal fluctuations of a scaled Hamiltonian 
H = vanish with H — > oo. In algebraic systems, 

this parameter needs to be introduced on the level of in- 
dividual generators, via scaled generators Gi = Gi/H K 
(with k > 0) whose substitution into the Hamiltonian 
H{Gi) should be consistent with the definition of H, 
thus H(G l /^ K ) = H(Gi)/X. In fact, this is why the 
size parameter H = was included into the effective 

coupling constant X/V^ of the above Hamiltonians 
(n = 1,2,3). The known commutation relations for the 
bare generators Gi then ensure that the scaled genera- 
tors Qi yield vanishing commutators in the H — > oo limit, 
[Gi,Gj] — > 0, which constitutes the classical behavior of 
the correctly scaled observables. 

A general method for approaching the classical limit 
in finite algebraic models is based on coherent states 
[33l Hi|. These for the above-described composite sys- 
tems are naturally considered in the form of a tensor 

product |C) ® 10) where 10 «■ e Cb+ |0) (with ( e C) is 
the HW(1) coherent state of the subsystem (i), and |0 is 
a yet unspecified SU(1,1) or SU(2) coherent state of the 
subsystem (ii) (3|1 Hg|. The latter states can be taken 
in several alternative forms, depending on a concrete re- 
alization of the two algebras. One possibility is to use 
|0 oc e iK +\k,0) or |0 oc e^ J +\j,-j) (with f e C) and 
associate the classical limit with j — > 00 or k 00. 
The corresponding phase space of subsystem (ii) is then 
identified with a 2D surface of constant positive or neg- 
ative curvature, which is the sphere j 2 + j 2 + y\ = 1 in 
the SU(2) case (ji = J%/j) or a two-sheet hyperboloid 
kl + k 2 ,- k 2 = -1 in the SU(1,1) case (fc 4 = K t /k). 

We do not directly follow this path, partly because in 
our case the value of the SU(1,1) invariant is fixed to 
k = \ or |, so we would not be able to keep the same 
treatment in both SU(1,1) and SU(2) systems. Instead, 
we employ the Holstein-Primakoff transformation of both 
algebras onto a bosonic field c 1 , c, which reads as 

K+ = J(2k + Jc)*, K_ = (2fc + c t c)5 c , 

K = Jc + k, (12) 

in the SU(1,1) case (with c f c = N c > 0) and 

J+ = c){2j - ch)i, J- = (2j - c)c)^c, 

Jo = c*c-j, (13) 

in the SU(2) case (with < N c < 2j). In this way, we 
obtain a mapping of the original HW(1) <E) SU(») dynam- 
ical algebra onto a new algebra associated with both b 
and c types of bosons, which is then analyzed with the 
aid of bosonic coherent states |C,0 ^ e |0), where 

|0) is a common vacuum of both b and c bosons. 



For Hamiltonians ([6]) and ([9]), the new dynamical 
algebra can be identified with the algebra U(2) = 
{b*b, e'e, b^c, c^6}, since in both these cases the total 
number of bosons N = Nb + N c is conserved. The size pa- 
rameters AfW and introduced above coincide both 
with the value 2N. On the other hand, for the noninte- 
grable Hamiltonian (|10l) the new dynamical algebra can 
be identified with HW(2), the Heisenberg-Weyl algebra 
of b and c bosons. The total number of bosons is not 
conserved and therefore the only sensible size parameter 
is the value of j (for consistency reasons we have chosen 
M< 3 ) = 4j), which measures the size of the subsystem 
(ii). 

The method proceeds via evaluating the expectation 

value (C,£Lff (n) |C,0/ M(n) for an y 01 the above Hamil- 
tonians H( n > (with n = 1,2, 3) in the coherent states of 
b and c bosons. The result can be gained directly by 
substituting 

c _ x + ip b _ y + iq 

and the Hermitian conjugate expressions for and 
into the scaled Hamiltonian U {n) = H^/M^ . Note 
that in Eq. (fT4]l we define coordinates x, y and the asso- 
ciated momenta p, q, respectively, which satisfy canoni- 
cal commutation relations [x,p] = [y,q] = i/M^ n ' with 
l/M^™) playing the role of Planck constant. 

For M( n ) — ¥ 00, the coordinate and momentum oper- 
ators can be treated as commuting variables. We obtain 
a general form 

H (n) = n (n) + m ,(n) ^ 

for the scaled classical Hamiltonian. The first term, 

nj(n) R^LOo UJq 2 2 U) 2 2 

^0 = 2 — ' 2 ^ y ^ ' 

which describes the system with A = 0, is common to 
all three models n = 1,2,3, with the additive constant 
expressed as i?W = 2k/M ( - 1 \ = 2j/M { - 2 \ and 

R^ = 2j/M^ = 0.5 (cf. Tab.Hl. The second term cor- 
responds to the interaction and has a mo del- dependent 
form 

= yj 2i?« + (p 2 + x 2 ) (xy + pq)/V2 , 

n' {2) = y/ 2i?( 2 ) - (p 2 + x 2 ) {xy +pq)/V2, (17) 

H' {3) = \/2R& - (p 2 +x 2 )V2xy. 

The constraints (J?} and (TTT|) on the conservation of 
M (1) and Af (2) read as 

p 2 + q 2 + x 2 + y 2 = 1 , (18) 

which makes it possible to completely eliminate one de- 
gree of freedom in both integrable models. To do this, 
we set one of the momenta to zero, in our case q = 0, 
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and use Eq. JT8 



to fix the corresponding coordinate: 
y = ±yl — p 2 — x 2 . On the level of coherent states, 
the choice q = is achieved by considering only a rel- 
ative phase between b and c bosons, setting the overall 
phase factor to unity. Note that this choice is dynam- 
ically consistent since the elimination of y ensures that 
q = dU/dy = 0. 



B. Ground-state phase transitions 

To analyze the H — > oo properties of the ground state 
of the above three models, we set both momenta p and 
q in Eqs. (fTo) and (IT71) to zero, yielding a potential 
y(n) = H {n) \p= q=0 (n = 1,2,3). Indeed, any increase of 
p or q takes one away from the minimum of the Hamilto- 
nian (p, q, x, y) and therefore corresponds to an exci- 
tation of the system above the ground state. The energy 
obtained by the minimization of the potential V^ n \x,y) 
represents an estimate of the scaled ground-state energy 
So = Eo/H. 

The problem can be further simplified for both inte- 
grable systems, i.e., for Hamiltonians ifW and H^ 2 \ 
where the constraint (fT8|) with p = q = restricts the 
ground-state solution to the unit circle x 2 + y 2 = 1. In 
the SU(1,1) case (n — 1) with uiq > ui, we take 



-0.5 



sin ■& , y — cos •& . 



(19) 



yielding the potential 
V«=V« + ^sin 2 



sin(2i?) v / 2i?( 1 ) + sin 2 d. 

(20) 

For the integrable SU(2) Hamiltonian (n = 2) with ui > 
ujq, we redefine the angle z9 so that 



x = cos i? , y = sin & . 



(21) 



and 



V< 2 > = + 



Aw 



, sin 2 i? + -= 

2 V¥ 



— COS 2 

(22) 



lu — u\ > 0, while Vq 1 ' 1 = (oj + 



In both cases Aw 

R^u )/2 and V< 2) = (1 - R^)lo /2. 

For n — 3, the nonintegrable SU(2) Hamiltonian H^> 
of the Dicke model with luq = u>, the constraint (|18[) is 
not applicable. We take polar coordinates 



arriving at 



x = r cos i? , y = r sin $ , 



V (s) = V< 3) + ~ r 2 +^= r 3 sin(2i?) V2J?( 3 ) - cos 2 0, (24) 
V2 



(23) 



with V^ 3) = -R^lj/2. 

It is clear that the AfW oo limit of the SU(1,1) 
model with k = ^, | gives i?W — > 0. On the other hand, 



1.5 



0.5 




1.5 



0.5 



-1 -0.5 

& [rad] 

FIG. 2: The potential energy V (1) (#) and V (2) (tf) of the inte- 
grable SU(1,1) and SU(2) models as a function of a rescaled 
coupling parameter g = A/(\/2Acj). The thick curve demar- 
cates the trajectory of the potential minimum, the dashed line 
indicates the saddle point position. 



in the integrable version of the SU(2) model it is natural 
to take M ( V = 4j, that is i? (2) = 0.5. Indeed, with this 
choice the total number of bosons N can be arbitrarily 
partitioned into Nj, and N c , including the extremal choice 
(Nb, N c ) = (0, N), which corresponds to the photon vac- 
uum combined with a fully excited array of atoms. With 
these settings, both expressions (f2"0")) and (|2"12j) become 
identical (except the additive constants). The potential 
energy surface alias V*- 2 -* is shown in Fig. [2] 

The minimum of both potentials V (1) and V (2) for A = 
is at $ — i?o(0) = (we may equivalently choose i?o = i"j 
which would have no influence on the conclusions below) . 
For increasing A, the minimum i?o(<^) remains at the same 
place until the critical value 



A 



(i) 



A 



(2) 



Auj 
V2 



(25) 



is reached. Here, •& — becomes a saddle point and the 
minimum i?q(A) deviates to negative values, following a 
trajectory 



sin 2 ti = 



12g 2 - 1 - V /I2g 2 + 1 _ A 



I8g 2 



1 

V2 Auj ~ 2 



(26) 

This nonanalytic evolution represents a second-order 
quantum phase transition. From Eq. (|26[) for g > g C Q — 
0.5 we get i?o(<?) ~ Vd — ScO which means that the criti- 
cal exponent for the order parameter i9q ( or = sin^o) 
is equal to \. 
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1 -1 



FIG. 3: The potential energy surface V^(x,y) of the Dicke 
model with u)o — uj. Values of a rescaled coupling strength 
g = X/{V2u) are given in each panel. 



For the Dicke Hamiltonian the dimension of the 
system is not reduced, so the properties of potential (|24|) 
must be analyzed in the plane (r, i?) = (x,y). Neverthe- 
less, the physics is similar to that described above. For 
A growing from to a critical value 



v(3) 
K c0 



ujujo 



(27) 



(A 



(3) _ 



j/V2 for ujq — ui) the potential minimum is lo- 



cated at (x, y) = (0,0), which corresponds to a separable 
state of unexcited atoms and the field vacuum. At the 
critical point (|27[) . the determinant of the Hessian matrix 
(composed from second derivatives of V' 3 ) with respect to 
both variables) evaluated in the minimum becomes neg- 
ative, which means that (x, y) — (0, 0) becomes a saddle 
point of the potential. Starting at this point, two de- 
generate minima deviate symmetrically to the quadrants 
with xy < (for A > 0). These minima correspond to 
a (nearly) degenerate parity doublet of the ground state 
solutions involving excitations of both atomic and field 
subsystems. The distance ro of the minima from the ori- 
gin increases with g = X/y/2uJoJo as y/g — g c o above the 
critical point g c o = 0.5, so we have again a second-order 
QPT with the order parameter ro characterized by the 
critical exponent i. Various stages of evolution of the 
potential V^ 3 -* are shown in Fig. [3] 



C. Excited-state phase transitions 

Any excited-state phase transition can be recognized 
in the dependence of a quantum level density p(£, A) on 
the scaled energy £. At the ESQPT point £ C (A) this 
dependence shows a nonanalyticity whose type enables 
one to classify the critical behavior in agreement with 
the standard typology of thermal phase transitions. The 
nonanalyticity of p(£ , A) is reflected by specific discon- 
tinuous features in the flow of scaled energy levels £i(X) 
through the boundary £ C (X) [H,[l7j]. 

The semiclassical theory of the level density [37] leads 
to the decomposition 



p(£,X) = p(£,X) + p(£,X), 



(28) 



where p and p represent smooth and oscillatory compo- 
nents, respectively. The oscillatory component can be 
expressed as a sum over periodic orbits in the general 
form p = cos(h~ 1 Sk + 4>k), where Ak and <fik stand 

for an amplitude and a phase shift of the kth orbit con- 
tribution, while Sk = § p ■ dx represents the action over 
this orbit. In the limit fi — > this part leads to infinitely 
rapid oscillations which cancel out if the level density is 
integrated over an arbitrary narrow interval of energy. In 
this limit, only the smooth component of Eq. (|2"8"|) is rel- 
evant. It is expressed via orbits of zero length, yielding 
the formula 

p(£,X) = (2nh)- f J S(£ -U{p,x,X))dpdx, (29) 

v v ' 

n(s,\) 

where % is the classical Hamiltonian depending in general 
on /-dimensional vectors of coordinates x and momenta 
p. The quantity Q(£,X)d£ represents a 2/-dimensional 
volume of the available phase space for the interval of 
energy (£,£+ d£). 

As follows from these considerations, in the systems 
with synonymous thermodynamic (H — > oo) and classical 
(h — > 0) limits any kind of nonanalyticity in the energy 
dependence of the classical phase space volume generates 
an ESQPT on the quantum level. Such nonanalyticities 
may follow, e.g., from the presence of the Hamiltonian 
stationary points [l5T - fl9l l38M4f| . Because the quantum 
microcanonical entropy is proportional to a logarithm 
of p(£,X), the dependence of the level density on the 
scaled energy represents a key for the ESQPT classifica- 
tion (consistent with the classification of the correspond- 
ing thermal phase transition). For instance, a jump of 
fi(£,A) at £ C {X) corresponds to a first-order phase tran- 
sition. It is known that in systems with / = 1 there exist 
ESQPT effects even stronger than those of the first-order 
type. These effects are associated with an infinite peak of 
fi(£,A), which shows up the corresponding dependence 
p{£ , A). The origin of this behavior is often found in a lo- 
cal maximum of the one-dimensional potential at energy 
£ C {X). On the other hand, a softer type of nonanalytic- 
ity in Q(£ ,A), like a discontinuous or infinite derivative, 
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FIG. 4: Level dynamics for the SU(1,1) model (left) and for 
the SU(2) integrable model (right) with M (1) = M (2) = 100 
and Auj — 1. The scaled energies were obtained by an exact 
diagonalization. The ESQPT above the ground state critical 
point A c o = 0.707 is apparent in the bunching of levels around 
critical energies E&' = 0.5 and fi 2 ' = 0.25, respectively. 



causes a continuous phase transition. This typically hap- 
pens in the systems with more than one degrees of free- 
dom [ItJ ■ In case of a discontinuous (n — l)th derivative 
(n > 2) the transition is of nth order, and for a singu- 
lar (infinite) derivative the transition has no Ehrenfest 
classification. 

The ESQPT effects are present in the spectra of all 
three models described above. We saw in Sec. IIII Bl 
that the ground-state QPTs are located at the critical 
points A^ from Eqs. ([2"5]) and 1(57)). It turns out that 

(n) 

for A > X c0 the singularity propagates into the excited 
spectrum. Let us first consider the integrable Hamiltoni- 
ans H^ 1 ' and both corresponding to an effectively 

one-dimensional configuration space. When the critical 
point is reached in these systems, the global minimum 
of the potential or V^ 2 ' changes into a saddle point, 
which remains present for all values A > A^ or Aq 2 ? . The 
saddle point represents a singularity of f2(£,A), causing 
the strongest type of ESQPT characterized by the infinite 
peak in the semiclassical level density. To see this, recall 
that in systems with / = 1 the integral in Eq. ([29]) is 
equal to the period r of the single (uniquely determined) 
classical orbit at energy £, hence fi(£ , A) = t(£ , A). If 
£ coincides with the energy of the x — (•& = 0) saddle 
point, the period becomes infinite because x = p = 
is a stationary point (x = p = 0) of both Hamiltoni- 
ans VS^{p,x) and W 2 '(j),x). Let us note that the same 
type of ESQPT is observed in systems with one quantum 
degree of freedom showing a local maximum of the po- 
tential, for instance in the Lipkin model [HI, H]], and 
many others [l(| [0} ■ 

A local increase of the level density at the saddle-point 




FIG. 5: Expectation values of N a (left) and J z (right) for 
individual states across the spectrum at A = 1.5. The left 
and right panels, respectively, correspond to the SU(1,1) and 
the SU(2) integrable models with M (1) = M (2) = 2000. The 
ESQPT is indicated by needlelike singularities located at the 
critical energies. 



energy £c = Vq in both integrable models, i.e., the 
SU(1,1) and Jaynes-Cummings models (n = 1 and 2, 
respectively), is demonstrated in Fig. [4] The two panels 
capture the evolution of quantum spectra for both models 
with the interaction parameter A, showing clear indica- 
tions of the ground-state QPT and its extension into the 
ESQPT on the right-hand side of the critical point, which 
for Auj = 1 is at A$ = A^ } = 0.707 (see Tab. [J). The 
calculation was done in a finite-size case, but it shows 
well pronounced precursors of the phase transitional be- 
havior. 

Additional ESQPT signatures are depicted in Fig. [SJ 
which shows expectation values of operators proportional 
to K — N c + k and Jq = N c — j in individual ex- 
cited states as a function of scaled energy £ for a fixed 
value of A = 1.5. Specifically, we consider the operator 
Na/M^ = (2K - |)/MW for the SU(1,1) model and 
Jo/j for the Jaynes-Cummings model. Note that these 
two observables act as order parameters of the respective 
standard QPTs: their ground-state expectation values 
change from (N a ) — to (N a )o > 0, and from (Jq)q = j 
to (Jo)o < j, as A crosses the critical point A c o. 

In both panels of Fig. [5l the respective system is well 
above the QPT critical point. The energy dependence of 
the respective expectation value (N a )g and (Jo)e shows a 
cusplike shape with a singularity localized at the ESQPT 
energy f| n ' (cf. Fig. 0}. The two shapes are mutually 
reversed: while for the SU(1,1) model, the expectation 
value drops sharply to the lowest value at the critical 
energy, for the SU(2) model it has a needle-shaped max- 
imum. This is connected with a singular localization of 
the semiclassical wave function for £ = £c^ at the saddle 
point of the potential, i.e., at x = for the SU(1,1) and 
x = 1 for SU(2) model (in both cases i? = Q). This implies 
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(N c ) £c = for the SU(1,1) case (hence (N a ) £c = or 1 for 
even or odd systems, respectively) and (N c )s c = 
for the SU(2) case (so (J Q ) £a = |M (2) - j). An analo- 
gous effect (explained by infinite dwell times of a classi- 
cal particle at the stationary point) is known from one- 
dimensional systems with a local maximum of the po- 
tential 38]. The ESQPT critical energies £ c (1) and 
drop to the ground-state energy as A decreases to the 
respective critical points xJq and A^q , and so do both 
cusp singularities in Fig. [5] Below the critical point the 
singularities disappear. 

For the nonintegrable Dicke model with the Hamilto- 
nian H^ 3 > , the two-dimensional potential (|24[1 has a sad- 
dle point at (,x, y) = (0, 0). This is connected with a non- 
analytic dependence of the phase space volume (and the 
level density) on the scaled energy, although of a softer 
type than in the previous case, as follows from a higher 
dimensionality of the phase space for the Dicke model. 
Specifically, for A > X y c0 ' the level density exhibits an 
anomalous growth with an infinite derivative (singular 

tangent) at £ = £ c , which coincides with the saddle- 
point energy of the potential [l7|. The resulting 
ESQPT is continuous (but without the Ehrcnfcst clas- 
sification), although its finite-size precursors very much 
resemble those of a first-order phase transition (the level 
density is close to a step-like function). 



The step-like increase of the level density in the Dicke 
model can be seen in Fig. [BJ where the level dynamics 
with variable A is shown for j — 2. Even for such a 
moderate value of the angular momentum, a sharp pre- 

(3) 

cursor of the ESQPT effect at absolute energy Ec = 
— f is well visible in the spectrum above 
A^o = 0.707 (for u) — uj = 1) as the lower interface 
between the horizontal and sloped level contours. Note 
that the effects of the Hilbert space truncation (the cutoff 
for the number of photons; see Sec. Ill Dj ) become relevant 
for the high-energy part of the spectrum. 



IV. QUENCH DYNAMICS 

A. Survival probability and energy distribution 

The Hamiltonians introduced in Sec. [TT1 have the com- 
mon form (|15p. that is H(X) = Hq + XH' if omitting the 
model specifying superscript n. Here "Ho an d W repre- 
sent the free and interaction term, respectively, and A is a 
dimensionless control parameter. Let us stress that here 
we are working with the scaled Hamiltonian H = -ff/K, 
but consider a finite-H case, so that in general [Ho, 'H'} ^ 
0. As seen from the expression "H(A 2 ) = 'H(Xi) + A'H' 
with A = A2 — Ai, the above Hamiltonian allows one to 
apply perturbation techniques with the same perturba- 
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tion H' for all initial points Ai. 

Suppose that the system is initially prepared in one 
of the eigenstates |^,-(Ai)) = of %{Xi) = Hi with 
energy i?i(Ai)/H = £\. Below we will consider the initial 
state l^i) coinciding with the ground state |^>o(^i))> but 
the formalism can be very easily developed for the general 
case. At time t = 0, the value of the control parameter 
is abruptly changed from Ai to A 2 = Ai + A. The state 
is no more an eigenstate of the new Hamiltonian 
H(X 2 ) = %2 and starts evolving. 

The evolution after the quench can be monitored by a 
survival probability pi(t) — |ai(i)| 2 , where 



ai(i) = (^i|e- m2t |Vi) 



\{£ 2 \i>i)\ 2 e- l£ * t d£ 2 

wi(£ 2 ) 



(30) 



is an amplitude describing the decay and recurrence of 
the initial state \ip\) for t > 0. A formula of this form 
captures in general all quantum decay processes and has 
been studied in many different contexts (e.g., in analyses 
of the fidelity or Loschmidt echo [Ilj]). Note that the use 
of the scaled Hamiltonian W 2 in Eq. ([30)1 is equivalent to 
the t — > f/N transformation of time in the expression with 
unsealed Hamiltonian H 2 . Expanding the initial state 
j^i) in the eigenbasis \£ 2 i) = \£i(X 2 )) of the Hamiltonian 
H 2 (with i = 1,2,... enumerating discrete eigenvalues 



|Vl>=^^2i|Vl>|f2i) 



(31) 



the survival probability reads as 
Pi(t) = £ N 4 + 2]T \ Cl \ 2 \ Cj \ 2 cos[(£ 2l - £ 2j )t] . (32) 

i i>j 

As indicated in Eq. (|30[) . the survival amplitude a\(t) 
can be written as the Fourier transform of the energy 
distribution wi{£ 2 ) = | (£2^1) 1 2 01 the initial state in the 
eigenbasis oiT-L 2 . The precise energy distribution is given 
by 



(£ 2 ) = J2h\ 2 ^2-£ 2i ) 



(33) 



Since both functions pi(t) in Eq. (|3"2")l and uji(£ 2 ) in 
Eq. ([53")) are expressed in terms of the discrete energies 
£ 2 i and the corresponding occupation probabilities |cj| 2 , 
they comprise fully equivalent information on the quench- 
induced relaxation process. 

The discrete form (|33p of the energy distribution 
u\{£ 2 ) can be approximated by its smoothened form 
lui{£ 2 ), obtained by replacing the 5 functions by normal- 
ized Gaussian profiles centered at eigenenergies £ 2 i. This 
leads to 



wi(£ 2 ) 



1 



: exp 



{£2 - £ 2i f 



2ai 



(34) 



where the widths <Ji are chosen separately for each Gaus- 
sian with regard to the local density of states in the re- 
spective part of the spectrum of T~i 2 . The aim of the 
smoothening procedure is to overcome the discrete char- 
acter of u>i(£ 2 ) while loosing as little as possible informa- 
tion on its local behavior. We therefore set the width of 
each Gaussian to the spacing between the («+l)th and 
ith levels, so <Ji = £2(2+1) ~ &2i- 

Although the discrete and smoothed forms (|33[) and 
(f34| capture basically the same information, their dis- 
tinction leads to two visualization methods of the en- 
ergy distribution. In the first one, based on the discrete 
form u)x{£ 2 ), the values of \a\ 2 are drawn against £ 2i in 
the form of a scatter plot (individual points being enu- 
merated by the eigenvalue index i). The other method 
shows the smoothened distribution uji{£ 2 ) as a continu- 
ous function of energy £ 2 . While the first method dis- 
plays essentially the average of \ci\ 2 in the given energy 
domain irrespective of the number of states (level den- 
sity) in this domain, the second method inherently con- 
tains a density-dependent weighting. In the following, 
we use both methods and compare the resulting forms 
with the time evolution of the survival probability pi(t), 
calculated from the exact formula (l32l). 



B. A critical quench 

Rather specific shapes of the energy distributions 
uii(£ 2 ) and u>i(£ 2 ) can be expected if the system exhibits 
an excited-state quantum phase transition for H — > 00. 
Assume that a sequence of such transitions is indeed 
present at energies £ C (X) depending, in general, on the 
control parameter A (we know that in the models stud- 
ied here, £ c is a constant). The critical curve £ C (A) in the 
plane £ x A may eventually reach the lowest energy of the 
system; then the ground-state quantum phase transition 
is observed at the corresponding value A = A c o (it is so 
in the present models). As discussed above, the flow of 
energy levels as a function of A and the density of the 
spectrum as a function of £ are nonanalytic when cross- 
ing the ESQPT critical curve. Therefore, if the energy 
distribution 0^1(^2) or bj\{£ 2 ) of the initial state after the 
quantum quench interferes with the critical value £ C (A 2 ), 
one may expect some anomalous properties of the sur- 
vival probability. 

We can easily estimate which parameter changes A = 
A2 — Ai may lead to such anomalous relaxation processes. 
To do so, recall that the mean value £ 2 of both energy dis- 
tributions uji{£ 2 ) and uji{£ 2 ) (both forms yield the same 
value) is given by £ 2 — (^llT^l^i), so 



£ 2 =£ 1 + &{i> 1 \'U'\i> 1 



(35) 



The "critical quench" A c for a given initial state is the 
one for which the average £ 2 coincides with the critical 
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FIG. 7: Graphical determination of the critical quench Ai — > 
A2 for a given initial state: the value of A c is given by an 
intersection of the tangent £[ with the critical "curve" £ c - 



value £0(^2), hence 



g e (Ai + A c ) - £1 
f 



(36) 



Of course, an actual range of the ESQPT-influenced 
quenches covers a wider interval of the A values around 
A c , depending on the width of the distribution 101(82) 
or uj\(£2) and also on the smearing effects in the phase- 
transitional signatures due to the actual finite value of 
the size parameter N. Nevertheless, the above formula 
yields a good estimate of a central point of the inter- 
val, where the quench dynamics can be expected to show 
strong ESQPT precursors. 

The application of the Hellman-Feynman theorem to 
Eq. (35]) yields (V'iI'H'I^i) = £'1 = d£ i {X)/dX\ x =x 1 ■ This 
leads to a simple graphical interpretation of Eq. (|36j) 
shown in Fig. [7J The parameter change Ai — > A2 is iden- 
tified with the critical quench if the tangent of the initial 
energy level £j(A) at A = Ai crosses the critical curve 
£ C (A) at A = A2. Let us stress that the final value of the 
control parameter corresponding to the critical quench, 
A2 = Ai + A c , differs in general from the critical value 
A c o of the ground-state quantum phase transition. 



C. Results for integrable models 

Now we are ready to discuss model-specific results for 
the energy distributions ([33]) and (|34| and the corre- 
sponding survival probability based on Eq. (13T)1) . Note 
that in this and in the following section we return to the 
unsealed energy E — M"?, 

We start with the two integrable models. Figure [8] 
shows results for three quenches in the SU(1,1) model 
with 2N b + N a = 2000 (thus N a even). The initial state 
is identified with the ground state at A — Ai = 1.5 and 
the respective final parameter value A2 is written sep- 
arately in each panel. In the upper row of panels we 



present the quantity u)\(E2) as a scatter plot of the val- 
ues \ci\ 2 versus the energy eigenvalue E 2 i- Note that the 
number of points is so large here that the scatter plots 
look like continuous curves. The panels from left to right 
correspond to a quench above, at, and below the criti- 
cal energy, which for the present setting coincides with 
E c = 1000. While for both noncritical quenches (left and 
right panels) the distribution of \ci\ 2 exhibits just a single 
peak centered at energy E 2 depending on the value of A2 , 
the critical quench to the final value A 2 = 0.936 (middle 
panel) leads to a more complex distribution. In this case 
we observe a double peak structure in the plot of |ci| 2 , 
the peak-separating minimum being localized exactly at 
the ESQPT energy. 

The criticality of the quench to A2 = 0.936 can also 
be seen in the other rows of panels in Fig. [5] In the sec- 
ond row, we show the smoothened distribution QxiEz) 
from Eq. . A clear difference from the first row is ob- 
served in the middle panel, where the second row shows 
a sharp maximum at the critical energy, in contrast the 
minimum in the first row. This is due to the above- 
discussed (Sec. IIV A|) distinction between the visualiza- 
tion methods based on the raw and smoothened energy 
distributions lji(E 2 ) and ui^E?). We know (Sec. MI C|) 
that there is a local increase of the level density around 
the critical energy E « E c (see Fig. @| connected with 
diverging periods of the classical trajectories passing the 
saddle point. This leads to a sizable increase of the distri- 
bution cDi(£ l 2), despite the fact that individual values of 
\ci\ 2 are lower in the critical region, as seen in the upper 
panel of Fig. [5] 

In the lower row of panels in Fig. [8] the survival prob- 
ability pi (t) is shown as a function of time for the three 
quenches discussed above. Again, similar patterns are 
observed for both noncritical quenches (left and right 
panels). In these cases, the survival probability exhibits 
regular damped oscillations. The time constant r of the 
decaying envelope is related to the total width AE of 
the associated peak in the energy distribution by the 
Heisenberg-like relation r oc 1/AE, while the frequency 
and form of particular oscillations depend on the mean 
energy and the fine structure of the energy distribution. 
For the critical quench (middle panel), the survival prob- 
ability behaves differently than for the noncritical cases. 
The quick initial decay is followed just by small random 
oscillations in the region pi(t) « 0, avoiding the slowly 
damped recurrences present in the other panels. This 
type of dynamics is connected with the above-discussed 
modified form of the energy distribution shown in the 
upper panels of Fig. [3] 

It needs to be stressed that we are dealing here with 
a finite system whose behavior is unavoidably quasiperi- 
odic. Hence, strictly speaking, the lack of recurrences 
seen in the middle low panel of Fig. [5] can only be tem- 
poral, as follows from the exact formula for the survival 
probability in Eq. (|3"2"j) . It is known, however, that in 
realistic situations the quasiperiodicity of quantum evo- 
lution on the long time scales is beaten by decoherence 
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FIG. 8: Energy distributions and survival probabilities for three quantum quenches in the SU(1,1) model (M"' = 2000). 
The initial state is the ground state at Ai = 1.5 and the final parameter values A2 are given in each panel. The leftmost and 
rightmost panels in each row correspond to quenches above and below the critical energy, respectively, while the middle panel 
depicts a quench to the critical region. Upper row: the energy distribution of probabilities \ci\ 2 , see Eq. (|33p . Middle row: the 
smoothened energy distribution from Eq. (|34|l . Lower row: the survival probability from Eq. (|30[) . 



effects [42j . The short and medium time scales addressed 
in the present calculations are therefore most substantial 
from the practical viewpoint. We interpret the observed 
difference in the character of the quench-induced relax- 
ation process as an important dynamical consequence of 
the ESQPT. 

Rather similar results are obtained also for the SU(2)- 
based integrable model. In Fig. |H1 we show the same 
quantities as in the previous figure, but for the Jaynes- 
Cummings model and only for the critical quench. The 
calculation was done with j = 500 (M^ = 2000), so the 
ESQPT critical energy E c — 500. Again, the initial state 
coincides with the ground state at Ai = 1.5 and the fi- 
nal parameter value A2 = 0.936, for which the results are 
shown, corresponds to the critical case. The energy dis- 
tribution wi(£ , 2) (a scatter plot of \ci\ 2 values) is given 
in the upper right panel of Fig. [51 the smoothened en- 
ergy distribution Wi(£^) hi the upper left panel, and the 
survival probability p\ (t) in the lower panel. We observe 
essentially the same behavior as in the middle column 
of panels in Fig. [5] This is not surprising since the two 
models have a rather similar structure. 

In summary, it is clear that the presence of an excited- 
state phase transition in the spectrum of both integrable 
models has major impact on the quench-induced relax- 
ation processes. We observe that the survival probability 
quickly decays and shows no recurrences in the medium 



time scale for the QQs which lead the system to the ES- 
QPT critical energy. This behavior gives a strong sup- 
port to the conjecture proposed (in connection with the 
Lipkin model) in Ref. [22j |. 

D. Results for the nonintegrable model 

The nonintegrable Dicke model shows more complex 
behavior than the two integrable models discussed above, 
and is also more difficult from the numerical point of 
view. We show in Fig. [10] the results obtained for the 
critical quench in this model. The system is defined by 
j = 40 (i.e., it contains 80 atoms) and by the resonance 
condition loq = u = 1. In this case, the absolute energy 
corresponding to the ESQPT is E c — —40. As in the 
previous cases, the initial state before the quench is the 
ground state at Ai = 1.5, while the final parameter value 
corresponding to the critical quench is A2 = 1.02. The 
arrangement of Fig. [10] is the same as that of Fig. [9] 

As discussed in Sec. Ill Dt the infinite dimension of the 
Hilbert space of the Dicke model requires to pay an ap- 
propriate attention to the convergence issues. In our 
case, the stability of results against the truncation of the 
Hilbert space was checked by varying the cutoff param- 
eter iVt runc for the number of b bosons (photons) until 
the convergence was reached for the quantities consid- 
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FIG. 9: The same quantities as in Fig. [8] but for the critical 
quench in the Jaynes-Cummings model with = 2000. 



ered. In practice, one has to perform several runs of the 
computation with increasing value of iVt runc and com- 
pare the results obtained in each run. It needs to be 
stressed that an optimal value of the cutoff parameter 
(satisfying a plausibly defined convergence criterion) de- 
pends on the quantities considered and particularly on 
the relevant range of energy. The calculations presented 
in Fig. [TU] (with 80 atoms) were done including all states 
of the photon field up to Nt — Nt Tunc — 220. The dimen- 
sions connected with these high particle numbers are at 
the limit of our present computing capabilities. 

Due to the nonintegrability of the Dicke model, the 
behavior observed in Fig. [TU]is partly different from the 
behavior of the same quantities in Fig. |H1 In particular, 
the energy distribution of the \ci\ 2 probabilities and the 
smoothcncd distribution loi^E^) exhibit much stronger 
fluctuations than those of the integrable models. For 
this reason, an additional smoothing procedure (different 
from the one described in Sec. IIV A|) has been applied 
to the result of the calculation. The smoothed energy 
distributions are presented in the upper panels of Fig. 1101 
As we see in the upper left panel, the local maximum of 
U)\{E2) at the critical energy, clearly observed in both 
integrable models, is lost. However, the main feature 
of both energy distributions, which is their split form 
around the ESQPT critical energy, is well reproduced. 

Concerning the survival probability after the critical 
quench (the lower panel of Fig. HUl) . it behaves similarly as 
in the above integrable models. It should be noted, how- 
ever, that for the integrable models the survival probabil- 
ity yields its characteristic "critical shape" (with no re- 
currences) only in a very narrow interval around the value 
of A2 defining the critical quench. In contrast, the Dicke 
model yields a much wider interval of critical-like relax- 
ation responses. Only for A2 far away from the critical- 
quench value (smaller or larger), the \ci\ 2 and wi(i?2) 




FIG. 10: The same as in Fig. [§] but for the nonintegrable 
Dicke model with j = 40. The curve in the upper-right 
panel is cut on the low-energy side because of the additional 
smoothening procedure. 



distributions receive their typical single-peak shapes and 
the survival probability p\ (t) gets the corresponding form 
of damped oscillations. 

We may conclude that the ESQPT in the Dicke model 
affects the quench dynamics in a qualitatively similar way 
as in the integrable models, although its fingerprints in 
various QQ-rclated observables are fuzzier than those dis- 
cussed in Sec. IIV CI The observed differences are partly 
due to the fact that the ESQPT in the Dicke model is 
of a softer type than those in the SU(1,1) and Jaynes- 
Cummings models (see Sec. MI C[) . Another reason for 
the softening is the chaotic nature of dynamics in the 
Dicke model [3l| . which, generically, has a tendency to 
obscure the ESQPT signatures [l7j|. 



V. CONCLUSIONS 

We have studied the phase diagram and the nonequilib- 
rium dynamics of three models describing the interaction 
of a single-mode bosonic field with an algebraic subsys- 
tem based on either the SU(1,1), or the SU(2) algebras. 

The existence of an excited-state quantum phase tran- 
sition in both integrable SU(1,1) and SU(2)-based quan- 
tum models is revealed for finite systems as a local peak 
in the level density, which in the thermodynamic limit 
transforms into a singularity. In the nonintegrable SU(2)- 
based model, the ESQPT leads to a step-like increase of 
the level density, which limits to a dependence with an in- 
finite derivative. The ESQPT manifests itself also in the 
expectation values of quantum observables that depict 
singularities at the critical scaled energy. These signals 
of the presence of an ESQPT open the possibility of us- 
ing the concept of order parameter and to resort to the 
Landau theory to characterize and classify them. 
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We have investigated the consequences of an ESQPT 
on the relaxation dynamics after a quantum quench. 
Starting from an initial state that we choose as the 
ground state of the system for a specific value of the con- 
trol parameter, a sudden change of the control parameter 
is applied and the relaxation process is followed by solv- 
ing the time-dependent Schrodinger equation. This is 
done either exactly (for the integrable models), or in a 
truncated space (for the Dicke model) , where the conver- 
gence issues are taken into account. 

Various relevant magnitudes related to the relaxation 
process after the quench are studied. In particular, wc 
analyze the survival probability of the initial state after 
the quench, which is closely related to the energy distri- 
bution of the initial state in the eigenbasis of the new 
Hamiltonian. We see a dramatic effect in the survival 
probability for the critical quench that drives the system 
from the initial ground state to the critical energy do- 
main associated with the ESQPT at the new value of the 
control parameter. This effect is studied separately in 
the three models used. 

The two integrable systems subjected to normal (non- 
critical) quenches display the typical pattern of collapses 
and revivals with a smooth decaying envelope, as follows 
from the single-peak forms of the respective energy dis- 
tributions. In contrast, the critical quench produces a 
sudden destruction of the survival probability followed 
by small random oscillations. This specific response is 
connected to a more complex shape of the energy distri- 
bution, showing a kind of splitting right at the ESQPT 
energy. 



A similar phenomenon is also observed in the nonin- 
tegrable Dicke model. However, due to level repulsion 
the survival probability is reduced in amplitude and the 
critical region is much broader. Away from the critical 
region, the survival probability and the energy distribu- 
tion behave in a similar way as in the integrable models. 

In spite of the differences observed in integrable and 
nonintegrable cases, we believe that relaxation dynamics 
offers clear signals of an excited-state phase transitions 
in the Dicke type of models. The character of these tran- 
sitions deserves more studies. 
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